Brought to you by:

The following article is Open access

The Doppler Flip in HD 100546 as a Disk Eruption: The Elephant in the Room of Kinematic Protoplanet Searches

, , , , and

Published 2022 June 28 © 2022. The Author(s). Published by the American Astronomical Society.
, , Citation Simon Casassus et al 2022 ApJL 933 L4 DOI 10.3847/2041-8213/ac75e8

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

2041-8205/933/1/L4

Abstract

The interpretation of molecular-line data using hydrodynamical simulations of planet–disk interactions fosters new hope for the indirect detection of protoplanets. In a model-independent approach, embedded protoplanets should be found at the roots of abrupt Doppler flips in velocity centroid maps. However, the largest velocity perturbation known for an unwarped disk, in the disk of HD 100546, leads to a conspicuous Doppler flip that coincides with a thick dust ring, in contradiction with an interpretation in terms of a ≳1 Mjup body. Here we present new ALMA observations of the 12CO(2–1) kinematics in HD 100546, with a factor of 2 finer angular resolutions. We find that the disk rotation curve is consistent with a central mass 2.1 < M/M < 2.3 and that the blueshifted side of the Doppler flip is due to vertical motions, reminiscent of the disk wind proposed previously from blueshifted SO lines. We tentatively propose a qualitative interpretation in terms of a surface disturbance to the Keplerian flow, i.e., a disk eruption, driven by an embedded outflow launched by a ∼10 Mearth body. Another interpretation involves a disk-mass-loading hot spot at the convergence of an envelope accretion streamer.

Export citation and abstract BibTeX RIS

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

The first exoplanets were detected indirectly through the stellar reflex motion (Mayor & Queloz 1995). Likewise, the gas velocity field in protoplanetary disks bears the imprint of planet/disk interactions, and molecular-line observations have emerged as an alternative approach to the identification of protoplanets. The first observability predictions, by Perez et al. (2015), linked embedded planets to kinks or wiggle-shaped deviations from sub-Keplerian rotation in channel maps (sky images in position–velocity data cubes, hereafter Iv ( x )). These predictions appear to match the data in HD 163296 and in HD 97048, where the most conspicuous deviations are interpreted in terms of the location and mass of the perturber through comparison with hydro simulations (Pinte et al. 2018a, 2019).

New hopes for protoplanet detection are thus coming from disk kinematics, with the interpretation of molecular-line data using hydrodynamical simulations. However, the simulations required to link the channel-map data to embedded bodies tend to idealize the systems. The unicity of a given planetary system configuration remains to be explored. It turns out that wiggles or kinks are almost ubiquitous within a disk, and also very frequently observed, as they are picked up in at least nine disks (Pinte et al. 2020) from the DSHARP ALMA Large Program (Andrews et al. 2018), despite these observations not being designed to trace small velocity structures. It may be that the ubiquity of velocity deviations reflects the global impact of gravitational interactions with the embedded planetary systems or even large-scale perturbations in gravitationally unstable disks (Hall et al. 2020; Paneque-Carreño et al. 2021). Another uncertainty is that a fraction of these structures stems from synthesis imaging limitations or from structure in the underlying optical depth. In any case, although velocity structures are expected from planet/disk interactions, it is not possible to unequivocally infer the location of embedded planets directly from the observed kinks or wiggles.

An alternative to pick up and locate embedded protoplanets has been proposed based on the velocity centroid maps, v( x ) = 〈v Iv 〉 = ∫dv v Iv /∫Iv dv. After the subtraction of the axially symmetric background flow, the velocity reversal along the planetary wakes, right at the location of the planet, should be observable as a Doppler flip in molecular-line maps (Pérez et al. 2018; Casassus & Pérez 2019), i.e., as an abrupt and localized change of sign in the velocity centroid of molecular-line maps.

Very deep and long-baseline ALMA observations of HD 100546 (up to 12 km; Pérez et al. 2020) revealed the largest known velocity deviations seen in channel maps. After subtraction of the background axially symmetric flow, these deviations translate into a conspicuous Doppler flip, whose amplitude reaches one-third of the Keplerian velocity (Casassus & Pérez 2019). The root of the Doppler flip is constrained within ∼5 au, at a stellocentric radius Rf ∼ 28 ± 3 au. The amplitude of the flip requires a rather massive body, ≳ 5 Mjup for an interpretation in terms of an embedded planet. However, as noted by Casassus & Pérez (2019), its coincidence with the massive dust ring of HD 100546's disk is in contradiction with a massive protogiant, as theory predicts the clearing of a gap by the protoplanet. The only Doppler flip reported in the literature, corresponding to the largest velocity deviations known, is thus the elephant in the room of kinematic protoplanet searches. It is pointless to search for planets causing small velocity deviations, if the largest velocity deviation known remains unexplained.

In this Letter, we present new ALMA observations of HD 100546 in 12CO(2–1) and adjacent continuum, at unprecedented angular resolutions (Section 2). We analyze the line kinematics in terms of 3D velocity deviations from an axially symmetric flow (Section 3) and discuss these results in terms of the possible origins for the kinematic structures (Section 4) before summarizing our results (Section 5).

2. Observations

The data presented in this Letter correspond to ALMA observations of HD 100546, carried out in program 2018.1.01309.S (P.I. S. Pérez). The source was acquired in six executions blocks, for a total of 9.2 hr on source. The array was in configuration C-10, with baselines ranging from 83 to 15,238 m, and 41–45 active antennas. A log of observations is provided in Table 1. Quasars J1107–4449 and J1145–6954 were used for flux and phase calibration, respectively.

Table 1. Summary of ALMA Observations

DateN Ant.Execution BlockToSAvg. Elev.Mean PWVPhase rmsBaselineMRS
   (sec)(deg)(mm)(deg)(m)('')
2019-06-0445uid://A002/Xdd3de2/X22e546842.61.17.683-152380.4
2019-06-0445uid://A002/Xdd3de2/Xa95554741.51.18.083-152380.4
2019-06-0545uid://A002/Xdd3de2/X147a558436.01.19.583-152380.4
2019-06-0541uid://A002/Xdd4cf3/X282555042.30.98.8236-152380.3
2019-06-0641uid://A002/Xdd4cf3/Xe0c556137.90.99.0236-152380.3
2019-06-2145uid://A002/Xddf4b5/X5b5558441.81.312.283-161960.3

Note. Summary of the new ALMA observations presented in this work. The table shows the date of the observations, the total number of antennas in each execution, the total time on source (ToS), target average elevation, mean precipitable water vapor column (PWV) in the atmosphere, antenna-based phase noise, minimum and maximum baselines, and maximum recoverable scale (MRS) for the used array configuration.

Download table as:  ASCIITypeset image

The correlator was set up to provide four spectral windows (spws). The continuum was sampled with spw 0, 1, and 2, at 233.0 17.0 and 219.0 GHz, each covering 2 GHz and thus yielding a total bandwidth of 6 GHz. The 12CO(2–1) line, at a rest frequency of 230.538 GHz, was sampled in spw 3 with 158.74 m s−1 channels, which we resampled to the local standard of rest into 160 m s−1 channels. The data were calibrated by staff from the North America ALMA Regional Center.

In this report, we focus on the disk kinematics and continuum emission of the central regions of the HD 100546 disk, along its inner ring, which is only ∼ 0farcs2 in radius. We report exclusively on the new long-baseline data and do not concatenate with the more compact configuration previously presented in Pérez et al. (2020). Keplerian rotation spans 35 mas in the time interval of ∼2 yr between the acquisitions of the present data and that of the intermediate-length configuration. This arc is almost 3× the finest angular scales reported here (about 0farcs12 or even 0farcs09, as measured in the central source assuming that it is unresolved, see Figure 1, or as estimated with one-third of the natural-weights beam, see Cárcamo et al. 2018) and almost twice those from Pérez et al. (2020). The concatenated continuum image of the ring is thus smoothed to angular scales two to three times coarser. However, the concatenated 12CO(2–1) data is much more sensitive to extended emission. Its presentation is postponed to a forthcoming article on the large-scale disk kinematics.

Figure 1.

Figure 1. Continuum at 225 GHz. (a): Model image Im obtained with uvmem, with an effective beam 0farcs018 × 0farcs012 at 26° (as estimated with a elliptical Gaussian fit to the central source, which appears unresolved). (b) Restored image, with a Briggs robustness parameter of 0, corresponding to a beam of 0farcs024 × 0farcs016 at PA 31.1. The noise in the image is 9 μJy beam−1.

Standard image High-resolution image

The data were self-calibrated using the standard procedure. Each execution block (EB) was self-calibrated independently. We set up CASA task gaincal to average whole scans (option "solint" set to "inf"). The self-calibration process provided an improvement of a factor of 1.5–2 per EB in tclean images using Briggs weights with robustness parameter r = 0.0, with typical dynamic ranges (i.e., signal-to-noise ratio) of ∼175 to ∼200.

Following the self-calibration procedure, we performed the final synthesis imaging with the uvmem package (Casassus et al. 2006; Cárcamo et al. 2018), in the same way as in Casassus et al. (2021). In the present application, uvmem produces a positive-definite model image Im and corresponding model visibilities ${V}_{k}^{m}$ that fit the visibility data ${V}_{k}^{\circ }$ in a least-squares sense, i.e., by minimizing

Equation (1)

where ωk correspond to the visibility weights and where the sum runs over all visibility data (i.e., without gridding).

For the restoration of the continuum image, shown in Figure 1, we use Briggs weights with robustness parameter r = 0, resulting in a clean beam 0farcs024 × 0farcs016/31.1, where we give the beam major axis (bmaj), minor axis (bmin), and direction (bpa) in the format bmaj×bmin/bpa. This restored image can be compared to that obtained by standard packages such as tclean in casa.

The channel maps of the 12CO(2–1) line, shown in Appendix B, were restored with Briggs r = 1 for a beam 0farcs038 × 0farcs024/0.0. In the case of emission that is much more extended than the beam, the level of aliasing is reduced with masks that isolate the signal. Here we used the same strategy as in Casassus et al. (2021), with Keplerian masks obtained with the tool 12 developed by Teague (2020).

Moments maps were extracted with double-Gaussian fits to each line of sight, using the package GMoments 13 (Casassus et al. 2021). The line profile for a given line of sight is modeled as ${I}_{v}={\sum }_{j=1}^{2}{I}_{j}^{A}\exp \left(-\tfrac{{\left(v-{v}_{j}^{\circ }\right)}^{2}}{{\sigma }_{j}^{2}}\right)$. The top side, which faces the observer, is approximately traced with the brighter Gaussian, with the largest amplitude. The moment maps that result from the double-Gaussian fits are summarized in Figure 2. The following analysis of disk kinematics is based on Figure 2(c). Systemic velocity was set to 5.71 km s−1, based on the kinematic analysis presented below (Section 3).

Figure 2.

Figure 2. Moment maps in 12CO(2–1) using double-Gaussian fits in velocity. (a) The two-Gaussian velocity-integrated intensity. (b) The velocity dispersion extracted using the two Gaussians (this image has been median filtered). (c) The velocity centroid of the brighter Gaussian. (d) The amplitude of the brighter Gaussian, ${I}_{1}^{A}$. (e), (f) Same as (c) and (d), but for the fainter Gaussian (with amplitude ${I}_{2}^{A}$, these images have been median filtered). In (e), the continuum ring is overlaid against ${I}_{2}^{A}$ in a single contour at 1/20 peak and is roughly outlined as a region of fainter amplitude.

Standard image High-resolution image

3. Disk Rotation Curve and Kinematics

Disk rotation can be approximated with an axially symmetric velocity field $\tilde{{\boldsymbol{v}}}({\boldsymbol{r}})$ in terms of a 3D disk rotation curve $\tilde{{\boldsymbol{v}}}(R)$, with $\tilde{{\boldsymbol{v}}}({\boldsymbol{r}})=({\tilde{v}}_{R}(R),{\tilde{v}}_{\phi }(R),{\tilde{v}}_{z}(R))$ in cylindrical coordinates. The deviations from the axially symmetric flow are then ${\boldsymbol{u}}({\boldsymbol{r}})={\boldsymbol{v}}({\boldsymbol{r}})-\tilde{{\boldsymbol{v}}}({\boldsymbol{r}})$, where v ( r ) is the 3D velocity field of the gas. If the emission stems from a thin layer at the unit-opacity surface, the first moment of the axially symmetric flow along a line of sight $\hat{s}({\boldsymbol{x}})$ is

Equation (2)

where x = f PA,i,ψ (R, ϕ). The coordinate transform f depends on the disk position angle PA and inclination i and relates x to the point, seen in the disk's cylindrical system, at the intersection of the line of sight with the disk surface, which is approximated locally as a cone with opening angle ψ. The observed velocity deviations are thus $\delta {v}_{\circ }({\boldsymbol{x}})={v}_{\circ }-{v}_{\circ }^{m}$ and derive from the intrinsic velocity deviations u ( r ).

Here we apply the ConeRot package, 14 described in Casassus & Pérez (2019) and Casassus et al. (2021), to derive the 3D rotation curve $\tilde{{\boldsymbol{v}}}(R)=({\tilde{v}}_{R}(R),{\tilde{v}}_{\phi }(R),{\tilde{v}}_{z}(R))$, in disk-centered cylindrical coordinates. Details on the extraction of this 3D rotation curve are given in Appendix C. We focus on the vicinity of the Doppler flip and postpone an analysis of the whole disk kinematics to a forthcoming article. Although ConeRot can also estimate the disk orientation based on the line data alone, these orientation estimates are partly degenerate with the rotation curve. We therefore fixed the disk orientation to that obtained from the continuum, which was extracted using the MPolarMaps package (git@github.com:simoncasassus/MPolarMaps.git; Casassus et al. 2021). This package minimizes the variance in the radial profile of the continuum intensity to estimate PA and i. The result is PA = 323 ± 0fdg7, i = 41.7 ± 0fdg4. The ring center was estimated to be offset from the phase center (itself pointed at the star) by Δα =0farcs008 ± 0farcs0008 and Δδ = 0farcs001 ± 0farcs0007, and the inferred systemic velocity is ${v}_{{\rm{s}}{\rm{y}}{\rm{s}}}=5.71\pm 0.02\,{\rm{k}}{\rm{m}}\,{{\rm{s}}}^{-1}$

The rotation curve for the HD 100546 disk, with meridional flows, is shown in Figure 3. The rotation curve for a purely azimuthal flow is very similar to Figure 3(a), although with a somewhat higher stellar mass (2.13 ± 0.04 < M/M < 2.26 ±0.04). The range of possible values for the stellar mass, including both cases with and without meridional flows, is 2.12 ± 0.04 < M/M < 2.26 ± 0.04, where the upper limits stem from vertical Keplerian shear, while the lower limit stems from cylindrical rotation (see Appendix C). These mass estimates, which assume a distance of 108.1 ± 0.4 pc (Gaia Collaboration et al. 2018), bracket the photospheric estimate of ${M}_{\star }={2.18}_{-0.17}^{+0.02}\,{M}_{\odot }$ (Wichittanakom et al. 2020). The lower-mass estimate of 1.83 ±0.01 M reported in Casassus & Pérez (2019) was affected by the extraction of the moment map using a single Gaussian, which is more contaminated by the back side of the disk than the two Gaussians and which led to higher inclinations (of ∼ 46°; the stellar mass is sensitive to the ∼2nd power of inclination, Casassus et al. 2021, their Equation (10)), along with a missing $\sin (i)$ correction.

Figure 3.

Figure 3. Rotation curve in HD 100546. The regions in cyan correspond to the total extent of the bright continuum ring. From top to bottom, we show (1) the azimuthal rotation curve ${\tilde{v}}_{\phi }(R)$. The dashed horizontal lines are comparison Keplerian profiles with the corresponding stellar mass. The curve labeled "mid" is an extrapolation to the disk midplane assuming vertical Keplerian shear. (2) The vertical velocity component curve ${\tilde{v}}_{z}(R)$, where ${\tilde{v}}_{z}\gt 0$ points away from the disk midplane. (3) The radial velocity component ${\tilde{v}}_{r}(R)$, where ${\tilde{v}}_{r}\gt 0$ points away from the star. (4) The opening angle of the cone tracing the unit opacity surface for 12CO(2–1).

Standard image High-resolution image

Given the disk rotation curve, we now turn to the nonaxial kinematics, i.e., to local deviations from axially symmetric flow. Figure 4 presents a summary of the kinematics, as seen on the sky, and compares them with the continuum emission. The same kinematic structures are shown in a face-on deprojection in Figure 5. The deprojection technique accounts for the thickness of the CO(2–1) surface and is described in Casassus & Pérez (2019).

Figure 4.

Figure 4. Nonaxial kinematics in the disk of HD 100546. (a) The velocity centroid v from Figure 2(c). (b) v after subtraction of the axially symmetric flow, ${v}_{\circ }^{m}$, including radial and vertical components (i.e., meridional flows) in the rotation curve. (c) ${v}_{\circ }-{v}_{\circ }^{m}$, same as (b) but without meridional flows. The single red and blue contours correspond to ±0.7 km s−1. (d) Continuum from Figure 1(a), with an overlay of the contours for ${v}_{\circ }-{v}_{\circ }^{m}$ from (c).

Standard image High-resolution image
Figure 5.

Figure 5. Face-on view of the nonaxial kinematics in the disk of HD 100546. Annotations follow from Figure 4.

Standard image High-resolution image

4. Discussion

The velocity deviations from a purely azimuthal flow, inferred from Figures 4(c) and 5, are very similar to those reported in Casassus & Pérez (2019). There is a coherent velocity structure in the form of an arc that appears to coincide with the bright submillimeter continuum ring. The correspondence with structures in the continuum is shown in more detail in Figure 7, where we also see that the center of the flip coincides with a bifurcation along this ridge.

The Doppler flip extends over ∼120° or ∼56 au in azimuth, at a fixed orbital radius. This is best seen in the polar expansion of Figure 6, which also shows that the intrinsic velocity deviations u must reach ±1/3 Keplerian. As noted by Casassus & Pérez (2019), along the disk minor axis, the line of sight of the azimuthal components cancels, ${u}_{\phi }\hat{\phi }\cdot \hat{s}=0$. This implies that ur or uz contribute significantly to δ v( x ) in the blue side of the flip, at least along the disk minor axis. Likewise, along the disk major axis, the line-of-sight component of a radial component cancels, ${u}_{R}\hat{R}\cdot \hat{s}=0$, which implies that the red side of the flip is either due to uz < 0 (receding toward the back side of the disk), or to uϕ > 0, accelerating beyond Keplerian rotation at the ascending node.

Figure 6.

Figure 6. Top: polar expansion of δ v, from the face-on deprojection in Figure 5, and with a 1D rotation curve (no meridional flows). The horizontal dashed line corresponds to Rf = 28 au. The vertical yellow and green lines each indicate the directions of the disk major and minor axis, respectively. Middle: extraction of δ v at Rf = 28 au, δ v( x = f PA,i,ψ (Rf , ϕ)). Bottom: intrinsic velocity field u , assuming that the red side of the flip is entirely due to azimuthal velocity deviations, while the blue side is due to either radial or vertical components, such that uz > 0 corresponds to a wind, and uR < 0 to stellocentric accretion. The x-axis of all plots corresponds to offset azimuth and is oriented in the direction of disk rotation.

Standard image High-resolution image

With an extension of the rotation curve to 3D, including meridional flows, the blue side of the flip essentially disappears from the velocity deviation map of Figure 4(d), while the red side remains. The above qualitative interpretations are thus confirmed on average, not just along the disk axes.

Interestingly, in the face-on view of the nonaxial kinematics shown in Figure 5(d), the blue and red sides of the flip align closely with the brightest continuum ridge. The velocity deviations appear to follow the peak continuum emission. Given the hints for an optically thick 225 GHz continuum (see Pérez et al. 2020; along with the reduction of peak amplitude ${I}_{1}^{A}$ and the outline of the ring in absorption in ${I}_{2}^{A}$), the bright ridge is likely a temperature feature rather than a local increase in density. Thus, the strongest velocity deviations coincide with the hottest continuum emission.

In what follows we consider the origin for the observed kinematics, taking into account the constraints on the intrinsic velocity field u . We do not consider the possibility of planet–disk interactions driven by planets inside the cavity, because the largest kinematic signatures of such bodies is in their immediate vicinity. The possibility of stellar companions inside the cavity is ruled out by the sparse-aperture-masking data presented in Pérez et al. (2020).

4.1. A Protoplanet Outflow

The velocity fields traced by the rarer CO isotopologues do not exhibit such strong velocity deviations as in 12CO(2–1) (Pérez et al. 2020). The velocity disturbance thus stems from high above the midplane (at the 12CO unit-opacity surface, with an aspect ratio 2–5 times larger than that of the thermal scale height, h1 ∼ 2–5 hr). On the other hand, the blue side of the flip corresponds to strong vertical velocity components. Thus, the Doppler flip in the disk of HD 100546 could perhaps correspond to an outflow launched by a compact body at the location of the flip, which gains in velocity with distance from the source, as in stellar winds.

A wind in the disk of HD 100546, approximately stemming from the location of the flip, has in fact been proposed based on unresolved high-velocity SO rotational line emission (Booth et al. 2018). In addition to the very distinct kinematics of the SO lines versus CO(3–2), sulfur-bearing molecules are known tracers of shocks and outflows (Booth et al. 2018), which release the sulfur locked in dust grains. In this scenario, SO should be observed downstream of the outflow and should coincide with the blueshifted part of the CO(2–1) flip. The hotter continuum emission may perhaps correspond to the shocked material.

The outflowing velocities seen in SO could be rooted in a jet launched by planetary accretion through a circumplanetary disk, as predicted by MHD models (Machida et al. 2006; Gressel et al. 2013). For the flip to be consistent with the geometry of an outflow, it is necessary that the blue side of the flip corresponds to material whose bulk flow lies along the disk vertical axis. If the outflow was aligned with the disk midplane, given the phase of the flip, its effect would be a Doppler flip opposite in sign. The red side of the flip should correspond to material flowing away from the observer and therefore escaping from the back side of the disk.

But the outflow scenario faces problems, as the continuum ring appears to be quite optically thick. The 12CO(2–1) channel maps seem to outline the ring, and the back side in 12CO(2–1) cannot be seen under the ring (Pérez et al. 2020). In addition, the blue side curves along a constant radius, thus following the trajectory of bound gas in a circular orbit, rather than the bipolar geometry of an outflow.

Finally, existing theoretical studies on protoplanet MHD polar outflows have considered massive bodies (Machida et al. 2006; Gressel et al. 2013). As mentioned above, a giant planet should displace larger dust grains from its orbital track. Midplane outflows from smaller planets in the super-Earth-mass regime have been produced by non-MHD hydrodynamical simulations (Kuwahara et al. 2019). Perhaps such bodies could act like outflow sources if their equatorial plane is inclined relative to that of the circumstellar disk. The impact of a magnetic field in the super-Earth-mass regime regime remains to be investigated.

4.2. A Disk Eruption

Perhaps the Doppler flip could be more accurately described as a disk eruption, a surface disturbance driven by an embedded protoplanetary outflow. This outflow would not puncture the 12CO layer and would only push molecular material both vertically and in the direction of rotation. In this hypothesis, the planetary outflow does not escape the disk, and instead deposits momentum in disk material that is predominantly in Keplerian rotation, which thus curves around in its orbit in a fashion that is reminiscent of high-altitude volcanic plumes.

We now consider what kind of outflow would be needed to drive the observed intrinsic velocity deviations through such a disk eruption. We hypothesize that the observed δ v are compatible with the flow resulting from an expansion impulse Γ but that this flow does not puncture the surface. This impulse corresponds to the force per unit volume and mass caused by the ram pressure $\rho {v}_{s}^{2}$ inside an expanding shell with mass density ρ, width Δrs and radius rs ,

Equation (3)

A shell with mass Ms will not puncture the surface if its impulse is balanced by the external pressure $\rho {c}_{s}^{2}$,

Equation (4)

We avoid the question of the launching mechanism, and consider a shell blown at constant vs by an embedded wind or outflow, sustained with a mass-loss rate ${\dot{M}}_{s}$:

Equation (5)

If we evaluate ρ at the CO(2–1) unit-opacity surface H1 with a hydrostatic density profile,

Equation (6)

where f1 = H1/H ∼ 2−5, cs = ΩK H, and Σg is the gas surface density. For the outflow to have an appreciable impact, we require rs H. At the center of the flip in the disk of HD 100546 (M ∼ 2 M), we have a cylindrical radius Rf ∼ 28 au and a disk aspect ratio given by ConeRot of $h=\tan (\psi )=0.2$ for ψ = 11° (which gives a midplane temperature Tm = 44 K if f1 = 4, for a mean molecular weight μ = 2.17). If we set vs vf ∼ 1 km s−1, then $\tfrac{{\dot{M}}_{s}}{{{\rm{\Sigma }}}_{g}}\sim {10}^{14}\,$ cm2 s−1. In other words, for Σg ∼ 5 g cm−2, an outflow with ${\dot{M}}_{s}\sim {10}^{-6}\,{M}_{\oplus }\,{\mathrm{yr}}^{-1}$ should have an appreciable impact.

If we adopt a standard ejection to accretion mass ratio of ∼10 (Konigl & Pudritz 2000), we see that the growth of a super-Earth over 1 Myr could drive such a disk eruption at an altitude over the disk midplane of 4H. It is interesting to note that the above argument, if extended to f1 = H1/H = 5 (so with a midplane temperature of 28 K), yields a much smaller mass-loss rate of ${\dot{M}}_{s}\sim {10}^{-8}\,{M}_{\oplus }\,{\mathrm{yr}}^{-1}$, which would correspond to the accretion of a few percent of a super-Earth mass over 1 Myr. We conclude from this argument that while embedded outflows would have little impact in the denser midplane, they should lead to conspicuous velocity deviations at higher altitudes, which is consistent with the basic properties of the velocity deviations seen in the disk of HD 100546.

4.3. Disk Mass-Loading Hot Spot?

Could the observed nonaxial kinematics be caused by an accretion hot spot onto the surface of the disk at the location of the flip? Garufi et al. (2022) report localized SO and SO2 line emission right at the convergence of the proposed accretion streamers onto the disks in HL Tau and DG Tau, where they also observe strong non-Keplerian kinematical signatures. HD 100546 is a system with a known envelope (Grady et al. 2001) that has been linked to the idea of cloudlet capture as a secondary (or tertiary) ongoing accretion phase (Dullemond et al. 2019). Even though there is no relevant detection of accretion streamers in the CO(2–1) data around the region of interest, some features show similarity to the systems where such asymmetric infall has been proposed, as some spirals have been observed both in the optical (Grady et al. 2001; Ardila et al. 2007) and near-infrared (Avenhaus et al. 2014; Sissa et al.2018). On top of the spiral structure seen at large scales of several arcseconds, Avenhaus et al. (2014) additionally discovered a spiral in the inner system for which the base roughly coincides with the observed kinematic deviation presented here. Lesur et al. (2015) and Hennebelle et al. (2017) showed that asymmetric accretion from a surrounding envelope can produce unstable accretion shocks in the disk, from which spirals are generated and propagate through the disk. How the kinematics look close to the shock location still remains an unexplored topic of interest.

We discuss the possibility of the Doppler-flip signal being linked to such an accretion hot spot: In this scenario, the red part of the Doppler flip traces vertically infalling material that locally heats the disk and deposits its momentum. An accretion shock could be an explanation for the localized detection of SO emission (Booth et al. 2018), as thermal desorption produces hot SO gas molecules in the postshock area (Aota et al. 2015). In the scenario proposed here, the postshock area would coincide with the blueshifted part of the kinematic signal; the SO detection by Booth et al. (2018) is also strongly blueshifted. Whether a localized hot spot of infalling material can also be consistent with the blueshifted side of the signal and whether it can dynamically perturb the dust structures at the midplane is an important question that requires future numerical investigation.

5. Conclusions

In this letter we presented new observations of HD 100546 in 12CO(2–1), and adjacent continuum, at unprecedented angular resolutions, which we interpreted in terms of 3D velocity deviations around the axially symmetric flow driven by Keplerian rotation. Our main conclusions are:

  • 1.  
    The previously reported Doppler flip stands out as a large velocity deviation from a purely tangential and axially symmetric flow, of ±1.9 km s−1, corresponding to intrinsic velocity deviations of ±1/3 Keplerian (Figure 6)
  • 2.  
    However, the blue side of the Doppler flip disappears when accounting for vertical and radial flows in the axially symmetric background flow (Figures 4 or 5). Given the length of the blue side, which extends 90° in azimuth along a fairly constant radius, the bulk of the velocity deviations should thus be vertical (because the height of the emitting surface is fixed at the τ = 1 surface).
  • 3.  
    The nonaxial velocity deviations are very closely aligned with the brightest ridge along the continuum ring, and the center of the flip coincides with a bifurcation along this ridge (see Figure 7).
  • 4.  
    The central mass is 2.1 < M/M < 2.3, in agreement with photospheric mass estimates (Figure 3).

Figure 7.

Figure 7. (a) Sky view of the Doppler flip in the disk of HD 100546, as in Figure 4(d), but overlaid on a version of the uvmem model image that has been processed with an unsharp mask (as described in Pérez et al. 2020). The gray contours, at 6%, 7%, and 8% maximum, are chosen to highlight structure along the ring. (b) Same as (a), but from a face-on perspective.

Standard image High-resolution image

Pending dedicated hydrodynamical simulations, we tentatively propose a qualitative interpretation for the Doppler flip in terms of a surface disturbance to the Keplerian flow, i.e., a disk eruption, driven by an embedded outflow launched by a ∼10 MEarth body. Such an outflow would deposit vertical momentum in disk material that is predominantly in Keplerian rotation and would pollute the disk material in sulfur-bearing species downstream from the source of the outflow, thus accounting for the SO line data. Another possible interpretation could involve a disk accretion hot spot at the convergence of a disk-mass-loading streamer infalling from the envelope.

Whichever the nature of the strong velocity deviations in the disk of HD 100546, the present data and analysis show that protoplanet searches based on disk kinematics should consider alternative interpretations to the non-Keplerian structures observed in molecular-line channel maps. The largest velocity deviations can be vertical or radial and seen on the disk surface, rather than azimuthal and in the disk midplane, as would be expected from the gravitational influence of an embedded and massive planet on its immediate vicinity.

We thank the anonymous referees for comments that improved the manuscript, and Richard Teague for interesting comments and suggesting Figure 10 (in the Appendix) S.C., M.C., and P.W. acknowledge support from Agencia Nacional de Investigación y Desarrollo de Chile (ANID) given by FONDECYT Regular grants 1211496, ANID PFCHA/DOCTORADO BECAS CHILE/2018-72190574, ANID project Data Observatory Foundation DO210001, ALMA-ANID postdoctoral fellowship 31180050, and FONDECYT Postdoctorado grant 3220399. This work was funded by ANID—Millennium Science Initiative Program—Center Code NCN2021_080. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01309.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Facility: ALMA. -

Software: MPolarMaps (https://github.com/simoncasassus/MPolarMaps, Casassus et al. 2021), uvmem (https://github.com/miguelcarcamov/gpuvmem, Cárcamo et al. 2018), ConeRot (https://github.com/simoncasassus/ConeRot, updated version described in this article).

Appendix A: Channel Maps

The observed channel maps for the 12CO(2–1) data are shown in Figure 8.

Figure 8.

Figure 8. Channel maps from the 12CO(2–1) datacube. The corresponding line of sight vLSR, in km s−1, is annotated in black in each image. The continuum from Figure 1 is outlined in contours. The beam is 0farcs038 × 0farcs024/0.0, and the noise is 1.02 mJy beam−1.

Standard image High-resolution image

Appendix B: Separation of the Front CO Layer and Moment Map Extraction

The extraction of the disk rotation curve from the velocity centroid map rests on the assumption that it corresponds to the front CO layer, as we need to convert each line of sight to a unique location in disk-centric coordinates. As explained in Section 2, we approximately separate the two τ = 1 CO(2–1) layers with a double Gaussian, where the brighter Gaussian traces the front side. While this approximation breaks down in some lines of sight, such as where the back side of the disk is exposed on the near side of the disk minor axis, a double-Gaussian fit is in any case a better representation of the observed profile than that obtained with a single Gaussian. In addition, in lines of sights where the separation of the front and back CO layers is larger than the line width (see Figure 3 of Pinte et al. 2018b, for an illustration of the two CO layers), the line profile is double peaked, which for a single-Gaussian fit would bias the velocity centroid in between the two velocity components.

The separation of the two CO(2–1) layers results in qualitatively different structures for ${v}_{1}^{\circ }$ and ${v}_{2}^{\circ }$. The position of the peak velocity centroid of the brighter Gaussian (${v}_{1}^{\circ }$), in absolute value, drifts away from the disk major axis with distance from the star, roughly toward the northeast, as expected for a flared surface (see Figure 2(c)). The opposite holds for the fainter Gaussian (${v}_{2}^{\circ }$), whose peak velocity drifts toward the southwest (see Figure 2(f)). In addition, the continuum ring appears to coincide with a region of lower secondary Gaussian amplitude (${I}_{2}^{A}$; see Figure 2(e)) and can be outlined (roughly) in absorption. The separation of the two CO(2–1) layers, using the double-Gaussian fits, is illustrated in selected channels in Figure 9.

Example line of sights spectra are shown in Figure 10.

Figure 9.

Figure 9. Selected channel maps derived from the double-Gaussian fits to the 12CO(2–1) datacube of HD 100546 presented here. The label at the bottom right of each image correspond to vLSR in km s−1. From top to bottom, each row shows the observed maps (obs), then the double-Gaussian model image (mod), then that corresponding to the brighter Gaussian (g1), and finally to that of the fainter Gaussian (g2). The two Gaussian components show conspicuous differences, and each roughly corresponds to the front and back disk surfaces.

Standard image High-resolution image
Figure 10.

Figure 10. Example line of sight profiles in the HD100546 12CO(2-1) data for for two selected locations, with offset R.A. (α) and Dec. (δ) in arcsec as indicated. The locations for each line of sights are also plotted on the integrated line intensity map and on the second moment of the double-Gaussian fit.

Standard image High-resolution image

The continuum ring in the disk of HD 100546 has a rather small radius, of about 5 to 10 beams even in these long-baseline data, and the central cavity is devoid of CO(2–1) emission. In addition, the standalone long-baseline data of HD 100546 are not ideal to trace the larger scales. As a proof of concept for the use of the double-Gaussian fit to separate the two Gaussian layers, we have applied GMoments to archival data of HD 163296 in 12CO(2–1) (from the DSHARP ALMA Large Program (Isella et al. 2018), which we reimaged in the same way as for the HD 100546 data presented here). The intensity maps for each CO layer are shown for selected channels in Figure 11. The corresponding moment maps are shown in Figure 12, where the concentric continuum ring system of HD 163296 is seen in absorption in ${I}_{2}^{A}$.

Figure 11.

Figure 11. Selected channel maps obtained with double-Gaussian fits to the 12CO(2–1) datacube of HD 163296, reimaged from archival ALMA data. Annotations follow from Figure 9.

Standard image High-resolution image
Figure 12.

Figure 12. Moment maps in 12CO(2–1) toward HD 163296, obtained with double-Gaussian fits in velocity. Annotations follow from Figure 2.

Standard image High-resolution image

Appendix C: Extraction of 3D Rotation Curves with ConeRot

We use the same notation and reference frames as in Casassus & Pérez (2019). The sky frame is represented by ${ \mathcal S }$ and is oriented using Cartesian coordinates (x, y, z), with $\hat{y}$ aligned due north. A rotation about the vertical axis $\hat{z}$ into frame ${{ \mathcal S }}^{{\prime} }$ aligns $\hat{y}^{\prime} $ with the disk major axis, while the $(\hat{x}^{\prime} ,\hat{y}^{\prime} )$ plane is still parallel to the sky. A second rotation about axis ${\hat{y}}^{{\prime} }$, by an inclination angle i, matches the $(\hat{x}^{\prime} ,\hat{y}^{\prime} )$ plane with the disk midplane in frame ${{ \mathcal S }}^{{\prime\prime} }$.

If all of the emission originates from the top side of the disk, which is facing the observer, we have a bijection between the line of sight x and the polar coordinate of its intersection with the surface of the cone representing the disk surface, defined as the surface of unit opacity, with height H1(R) above the midplane. This cone has opening angle $\psi =\arctan ({h}_{1})$, where h1 = H1/R is the aspect ratio of the disk surface. We thus transform the sky coordinates x to cylindrical polar coordinates on the surface of a cone (R, ϕ), as measured in ${{ \mathcal S }}^{{\prime\prime} }$. The coordinate transform from Cartesian coordinates in ${{ \mathcal S }}^{{\prime\prime} }$ into cylindrical coordinates in ${{ \mathcal S }}^{{\prime\prime} }$,

Equation (C1)

depends on disk orientation and is invertible. Similar coordinate transforms have been used in, for example, Rosenfeld et al. (2013) or in Isella et al. (2018). We implement f with the following formulae, already provided in Casassus & Pérez (2019) but which we reproduce here for completeness:

Equation (C2)

Equation (C3)

The inversion of Equations (C2) and (C3) to obtain $(r,\phi )={{\boldsymbol{f}}}_{\mathrm{PA},i,+\psi }^{-1}({x}^{{\prime} },{y}^{{\prime} })$ can be achieved by noting that R is the root of

Equation (C4)

and with

Equation (C5)

If all of the observed signal stems from the top side of the disk and in the H1(R) surface, then the first moment of the axially symmetric flow along line of sight $\hat{s}({\boldsymbol{x}})$ is, using Equation (C1),

Equation (C6)

or

Equation (C7)

where all coordinates are measured in ${{ \mathcal S }}^{{\prime\prime} }$, i.e., in the frame of the disk. Similar formulae for the model velocity centroid have been proposed by Teague et al. (2019).

The origin of ϕ coincides with the PA of the disk major axis, defined as the line of ascending nodes (so the red side). The sign convention used here results in ${\tilde{v}}_{r}(R)\lt 0$ for accretion. In case of prograde rotation, with i < 90°, a wind flowing away from the midplane corresponds to ${\tilde{v}}_{z}(R)\gt 0$. But for retrograde rotation, the top side of the disk that is facing the observer corresponds to z < 0 in ${{ \mathcal S }}^{{\prime\prime} }$, so the disk opening angle ψ < 0, and for a wind ${\tilde{v}}_{z}(R)\lt 0$. The observed velocity deviations are thus $\delta {v}_{\circ }({\boldsymbol{x}})={v}_{\circ }-{v}_{\circ }^{m}$ and derive from the intrinsic velocity deviations u ( r ):

Equation (C8)

Given the disk orientation (PA, i, ψ), we solve for the disk rotation curve ${\tilde{{\boldsymbol{v}}}}_{\circ }(R)$ by performing a least-squares fit of ${v}_{\circ }^{m}(R,\phi )$ in Equation (C7) to v( f PA,i,ψ (R, ϕ)), which is the observed velocity centroid resampled by f PA,i,ψ . We obtain ${\tilde{v}}_{\phi }({R}_{l})$, ${\tilde{v}}_{r}({R}_{l})$, and ${\tilde{v}}_{z}({R}_{l})$ for each discretized value of Rl by minimizing

Equation (C9)

which assumes that the systemic velocity vs is known. In practice, we solve for the root of χ2 = 0 using the Cramer rule. The sum in Equation (C9) runs over all azimuths within a domain ${{ \mathcal D }}_{\phi }$, with weights

Equation (C10)

where Nb is the number of beam major axes along a circle with radius Rl , clipped so that Nb > 1. Nb approximately corrects for correlated pixels in. The domain in azimuth ${{ \mathcal D }}_{\phi }$ is typically [0, 2π], but it is interesting to restrict ${{ \mathcal D }}_{\phi }$ to [0, π], corresponding to the far side of the disk. The model velocity centroid ${v}_{\circ }^{m}(R,\phi )$ represents the axially symmetric flow and can be converted back to ${{ \mathcal S }}^{{\prime\prime} }$ by resampling with ${{\boldsymbol{f}}}_{\mathrm{PA},i,\psi }^{-1}$.

If vs is not known, then we initially restrict to a purely azimuthal rotation curve and fit for vs (Rl ) and ${\tilde{v}}_{\phi }({R}_{l})$ simultaneously, for all {Rl } in a relatively narrow radial domain [R1, R2], which we call a region. This optimization also includes the disk orientation (see next paragraph), but for constant orientation parameters (i, PA, ψ) representative of this region. The value of vs is then fixed to its median, and its uncertainty to its standard deviation, both extracted over overlapping regions that cover a wider range in radii. The full rotation curve is then calculated using this global vs value.

Disks are flared, i.e., the disk aspect ratio ${h}_{1}(R)\,=\tan (\psi (R))/R$ is a function of R. With molecular-line tracers, we expect that h1(R) will initially increase with R and then drop to zero as the tracers becomes optically thin at the outer edge of the disk. Disks may also warp, so that disk orientation will depend on R. We calculate the disk orientation profile (i(R), PA(R), ψ(R)) by performing the optimization in overlapping radial bins. The full radial extension of the disk is divided into M overlapping radial bins ${\{[{R}_{1j},{R}_{2j}]\}}_{j=1}^{M}$, thus defining radial regions $\{{{\rm{\Theta }}}_{j}({{\boldsymbol{x}}}^{{\prime} })\}{}_{j\,=\,1}^{M}$, where ${{\rm{\Theta }}}_{j}({{\boldsymbol{x}}}^{{\prime} })=1$ if $R={f}_{r}^{-1}({{\boldsymbol{x}}}^{{\prime} })\in [{r}_{1j},{r}_{2j}]$, and ${{\rm{\Theta }}}_{j}({{\boldsymbol{x}}}^{{\prime} })=0$ otherwise. In each radial bin, we maximize the log-likelihood function, −0.5χ2 to obtain (ij , PAj , ψj ) with

Equation (C11)

which is the radial sum of the weighted azimuthal variance of residuals relative to the expected variance ${\sigma }_{\mathrm{ref}}^{2}$, with σref typically of order 0.1 km s−1 when including image synthesis and channelization systematics. The sums extend over an interval in radius [R1, R2] and over all azimuths. We thus obtain an orientation profile ${\left\{i({R}_{j}),\mathrm{PA}({R}_{j}),\psi ({R}_{j})\right\}}_{j=1}^{M}$ as well as axially symmetric models $\{{v}_{j}^{m}({\boldsymbol{x}})\}{}_{j=1}^{M}$ in all radial regions. The axially symmetric model corresponding to these profiles is approximated by the average over all regions, i.e.,

Equation (C12)

This approximation ignores shadowing between regions along a line of sight, and is thus only applicable at low to moderate inclinations.

The minimization of Equation (C11) and the optimal orientation parameters, along with their associated errors, were achieved with a Markov Chain Monte Carlo ensemble sampler (Goodman & Weare 2010). We use the emcee package (Foreman-Mackey et al. 2013), with flat priors and with typically 300 iterations and 30 walkers.

The azimuthal rotation curve ${\tilde{v}}_{\phi }$ can be used to constrain the stellar mass under the assumption of pure Keplerian rotation. As noted in Casassus & Pérez (2019), the stellar mass is then bracketed by two extreme cases. Assuming a rigid vertical structure, in which the midplane velocity is the same as that at height H1, sets a lower limit. Pure vertical Keplerian shear, ignoring radial hydrostatic support, sets an upper limit, in which midplane azimuthal velocities are ${\tilde{v}}_{\phi }(R){\left(1+{h}_{1}^{2}\right)}^{3/4}$.

Footnotes

Please wait… references are loading.
10.3847/2041-8213/ac75e8